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ABSTRACT 

In order to assess the ability of helioseismology to probe the subsurface structure and magnetic field of sunspots, we need to determine 
how helioseismic travel times depend on perturbations to sunspot models. Here we numerically simulate the propagation of f, pi, and 
P2 wave packets through magnetic sunspot models. Among the models we considered, a ±50 km change in the height of the Wilson 
depression and a change in the subsurface magnetic field geometry can both be detected above the observational noise level. We also 
find that the travel-time shifts due to changes in a sunspot model must be modeled by computing the effects of changing the reference 
sunspot model, and not by computing the effects of changing the subsurface structure in the quiet-Sun model. For p t modes the latter 
is wrong by a factor of four. In conclusion, numerical modeling of MHD wave propagation is an essential tool for the interpretation 
of the effects of sunspots on seismic waveforms. 
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^1. A numerical approach 

$— i Elucidating the subsurface processes of sunspot formation 
^ and evolution is important for understanding solar m agnetism 
C3 ilMoradi et al.ll2010t iRempel & Schlichenmaier1l201 ll) . Helio- 
l— ~ 1 seismology is the only available tool for imaging th e subsurface 
I structure, magnetic field and dynamics of sunspots (iGizon et al.l 
^ 1201 Oh . all of which are largely unknown. 

i Early approaches interpreted the seismic signature of 
<^> sunspots in te rms of small perturbations to a qu iet-Sun ref- 
C~) erence model dFan et al 1 11995b lKosovichevlll996h . However, 
the effects of the surface magnetic field and associated struc- 
tural changes cause large perturbations in the wave speed (e.g. 
r^> iLindsev & Braunll2005l : IGizon et al.ll2006l) . Recently, direct nu- 
C^) merical simulations of MHD wave propagation through sunspot 
t-H models have been used to capture the large effects of the 
'<~J magn etic field on the waves (IGizon et al.1 12010): ICameron et al.1 

p\ . A study by iLiang et ail (120131) measured the travel-time 
H shifts caused by the sunspot in AR 9787. In an effort to un- 

- - -derstand and interpret these measurements we used numerical 
simulations. We simulated the interaction of surface-gravity (f) 
and acoustic (pi and P2) wave packets with several sunspot mod- 
els. With spectroscopy the height of the Wilson depression i n 
sunspots is not known to better than 50 km (ISolanki et al.ll 1993b . 
Here we investigate if seismic diagnostics can provide a stronger 
constraint. Further, we study the seismic consequences of two 
subsurface magnetic field geometries below sunspots: a mono- 
lithic flux tube versus a weaker, spread-out magnetic field. Fi- 
nally, we explore the sensitivity of waves to sound-speed pertur- 
bations at different depths along the axis of a sunspot model. 



These parametric studies enable us to test the assumptions 
of existing approximation methods used in the interpretation of 
sunspot seismology. By comparing with measured noise in seis- 
mic travel times we are able to draw conclusions about the de- 
tectability of different aspects of sunspot structure. 



2. Sunspot models 

2.1. Background quiet-Sun model (QS) 

Our sunspot models are embedded in t he background 
solar model CSM_A described by ISchunker et al.l 
d2011l). The CSM_A model is based on Model S from 



IChristensen-Dalsgaard et aD (Il996l) . but was modified to be 
convectively stable so that linear wav es can be propagate d 
using the SUM wave simulation code ( ICameron et al.l I2008I) . 
An indication of the sensitivity of the modes of oscillation in 
CSM_A as a function of height, z, is given by their kinetic 
energy density (see Fig. [TJ. Throughout this paper we define z 
as height measured upward from the bottom of the photosphere. 

2.2. Reference sunspot model (RS) 
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Ca meron et al 1 (l2011h constructed a parametric sunspot model 
that causes a seismic response consistent with that observe d for 
the sunspot in active region NOAA 9787 dLiang et al.l2013l , Pa- 
per I). The vector magnetic field of the sunspot model is cylindri- 
cally symmetric, i.e. it depends only on the horizontal distance 
from the axis of the sunspot, m — (x 2 + y 2 ) 1 ^ 2 , and on height, z. 
The model of ICameron et al.l d201 ll) is a self-similar and mono- 
lithic magnetic field model specified by the strength of the field 
on the axis of the sunspot (rn — 0) as a function of height. Here 
we adopt the same sunspot model with the parameters given in 
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Fig. 1. Reference sunspot model and perturbations to it. The mag- 
netic field lines of the reference model are shown by the black curves 
(Sect. 2.1.2). The dashed black curves are the magnetic field lines of the 
alternative magnetic field model (Sect. 2.1.4). Two Wilson depressions 
modified by ±50 km have been considered (red curves, Sect. 2.1.3). 
The ellipses along the axis of the sunspot indicate the locations of the 
various sound-speed perturbations (Sect. 2.1.5). The size of the ellipses 
corresponds to the full width at half maximum of the Gaussian pertur- 
bations. The right panel shows mode kinetic energy density, p||v|| 2 , for 
the f mode at kR Q = 450 (black), for at kR Q = 420 (blue), and for 
p 2 at kR e = 390 (red), where p is density and v is the velocity eigen- 
function for quiet-Sun model CSM_A. Mode kinetic energy density is 
a useful indicator for the vertical sensitivity of the modes. 



Table 1 of lCameron et al.ld201 lb . except that the parameters con- 
trolling the inclination of the magnetic field are modified slightly 
(ai = 1.5 Mm and = 9.2 Mm) so that the magnetic field incli- 
nation at the umbra-penumbra bound ary (50°) better match es ob- 
served field inclinations in sunspots ( ISchunker et al. 2008). The 
solid black curves in Fig.Q~|show a sample of the magnetic field 
lines. Throughout this paper we refer to this model as the refer- 
ence sunspot model abbreviated as 'RS'. 

2.3. Changes to the Wilson depression 

The Wilson depression is a physical depression of the surface 
(optical depth un ity) caused by strong mag netic fields in the um- 
bra of a sunspot (Loughhead & Bravll 19581) . The Wilson depres- 
sion of our reference sunspot model is at z = -550 km in the 
center of the u mbra and at z = -30 km in the middle of the 
penumbra as in lCameron et al. 

In order to assess the ability of helioseismology to measure 
the height of the Wilson depression, we consider two additional 
sunspot models in which the Wilson depression is shifted by 
+50 km with respect to the reference sunspot model. In these 
modified models, we translate the density, pressure, and sound- 
speed profiles in height by +50 km, while keeping the magnetic 
field unchanged (models are not in magnetohydrostatic equilib- 
rium). 



2.4. Alternative magnetic field model 

The monolithic reference sunspot model described in Sect. 12.21 
is only one possibility among many. We considered a sunspot 
model with an alternative magnetic field below z — -1 Mm, 



which fans out as a function of depth. This alternative model 
is motivated by, e.g . , the dynamical disconnection model of 
ISchlissler & Rempell (120051) whereby the flux tube is disrupted 
by convective motions at depths greater than 2 Mm. In the alter- 
native magnetic sunspot model, we keep the sound speed, pres- 
sure and density profiles unchanged from the reference sunspot 
model. 

More precisely, the alternative magnetic field is modified so 
that the magnetic field becomes weaker below z t = -1 Mm 
(above, it is identical to the reference sunspot model), smoothly 
transitioning to weak field strengths below Zb = -4 Mm. Along 
the axis of the sunspot, the z-component of the magnetic field, 
B z , is obtained from the reference field, B z , as follows: 



| 1 forz>z t , 
B z (z, vr = 0) = B z (z, m = 0)\ U(z) forz b < z < z t , 

I 0.1 forz <z b . 



(1) 



where II(z) is the fifth order polynomial in z such that B z (z) is 
continuous and has continuous first- and second-order deriva- 
tives, i.e. II(z) = 1 - 0.9Z 3 (6Z 2 - 15Z + 10) where Z = 
(z-Zt)/(zb-Zt). The same field lines as for the reference sunspot 
model are shown for the alternative magnetic field model in 
Fig. ID 

2.5. Subsurface sound-speed perturbations 

We also consider a set of sunspot models in which we modify 
the subsurface sound speed. In these models, the sound speed is 
given by 

c(tzr, z) — c(m, z) + Ac(w, z) (2) 
where c is the sound speed in the reference sunspot model and 

u 2 (z - hf 



Ac(vj, z) = c{vt, z)A exp 



2&L 



2cr? 



(3) 



is a 3D Gaussian perturbation to the sunspot model. In the above 
expression, the relative amplitude of the perturbation is A, the 
height of the centre of the perturbation is h, and the width of the 
perturbations are fixed to cr z = 0.42 Mm, cr w = 2.12 Mm. For 
one experiment (Sect. IBTTT i we use A — 0.1 and h = -2, -3 and 
-4 Mm. For a second experiment (Sect. l5T2t we use A = -0.01 
to 0.5 and h = -1.5 Mm. 



3. Numerical simulations of wave propagation 

We used the SUM code ( ICameron etal.ll2007l 1201 ll) to simu- 
late the propagation of wave packets through the sunspot mod- 
els described above. We use a simulation domain that covers 
145.77 Mm in each of the horizontal directions and extends from 
-25 Mm below to 2.5 Mm above the bottom of the photosphere. 
The depth of the simulation domain is sufficient to study acous- 
tic modes with radial orders « < 4. The spatial resolutions are 
0.025 Mm per pixel horizontally and 0.73 Mm per pixel verti- 
cally. The time resolution is 0.13 s. We did a resolution study 
which showed that for the waves studied here the numerical so- 
lutions are converged. 

We carry out simulations for individual wave packets with 
modes of fixed radial order, i.e. with either f, pi, or p2 modes. 
The initial conditions (f = 0) for the simulations are specified by 
the displacement, £(x,y,z), and the velocity, v(jc, y, z). This re- 
quires specifying the wave amplitude distribution, A^, as a func- 
tion of wavenumber, k (see lCameron et al .1200811201 ll) . For each 
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wave packet A* is chosen in such a way that the vertical veloc- 
ity in the simulation (v z ) is directly comparable to an observed 
cross-covariance function. This comparison is justified because 
the cross-covaria nce is closely related to the Green's function in 
the far field (e.g jGizon et"afll2010l) . 

The wave packets in our simulations are comparable to 
the Solar and Heliospheric Observatory's Michelson Doppler 
Imager ( SOHO/MDI) obser ved cross-covariance functions pre- 
sented bv lLiang et al.l d2013l) . The peak of wave power occurs at 
kR Q = 450 and v = 2.1 mHz for the f modes, at kR e = 420 
and v = 2.5 mHz for the pi modes, and at kR Q = 390 and 
v — 3.1 mHz for the p2 modes. At t = 0, each wave packet 
peaks at x — -43.73 Mm (at the left edge of Fig. [2J, from where 
they propagate in the +x direction (to the right) and interact with 
the sunspot centred at x — y — Mm. 

4. Travel times 

4. 1 . Measuring travel times 

Travel-time measurements are very useful to characterize the 
sensitivity of seismic waveforms to subsurface perturbations. 
In this paper we explore how changes to the reference sunspot 
model affect travel times. 

We define the travel-time shift, At, as the time by which a 
sliding reference wavelet, w, must be shifted to best match the 
seismic waveform v 7 ( x, y, f) - v?(x, y, z = 0.2 Mm, t). Following 
iGizon & Birchl (I2OO2I) and lLiangetai1(l2013l) At is the time, t, 
that maximises the function 

| v z (x, y, t) w(x, y, t — r) df 

F(t) = J -—^r- . (4) 

J \w(x,y,t)\ 2 dt 

In practice, the reference wavelet w was chosen to be v z from 
either the quiet-Sun simulation or the reference sunspot sim- 
ulation. To denote these two kinds of measurements (which 
give similar results, see Sect 15. 2t . we use the notations Aqs t 
and Arst respectively. Furthermore, we use the notation 
At(MODEL) to indicate the background model through which 
the seismic waveform is propagated, where 'MODEL' is, e.g., 
'RS', 'RS+pert', or 'QS+pert'. The time integrals in Eq. |4]are 
approximated by summations with a time discretisation of 6 s. 
With this time discretisation travel times are accurate to 0.01 s, 
which is sufficient for our purpose. 




-20 20 40 60 80 
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Fig. 2. Spatial map of travel-time shifts Aq S t for the reference sunspot 
model. The SLiM code numerically computes the waveform as a func- 
tion of time through the reference sunspot model. The color bar indi- 
cates the travel-time shift in seconds. The ellipse encloses the region for 
which the spatial average of the travel-time shift (denoted with symbols 
()) is computed. The grey (black) circle indicates the penumbra (umbra) 
of the sunspot. 

of the spatially averaged travel-time shifts. The observed noise 
level decreases with radial order, it is 3.8 s for the f modes, 3.5 s 
for pi modes and 1.6 s for P2 modes (Table[TJ. 



5. Results 

We carried out two sets of numerical experiments. In the first set, 
we propagated f, pi and P2 waves independently through differ- 
ent sunspot models to predict the travel-time shifts caused by 
various perturbations to the reference sunspot model, and to de- 
termine if they are above the observed noise level. In the second 
set of experiments we looked at the sensitivity of pi wave travel 
times to changes in the amplitude of the subsurface sound-speed. 
We considered sound-speed perturbations to both the quiet-Sun 
model and to the reference sunspot model. The results of these 
experiments are described in the sections below. 



4.2. Estimating travel-time noise 

While the numerical simulations performed in this paper are 
used to estimate the expectation value of a travel-time measure- 
ment (the signal), we also need to estimate the noise level in 
order to determine if a signal is detectable. 

We use the travel-time maps measured bv lLiang et al.l (1201 3l) 
covering se ven independent days of SOHO/MDI observations 
of AR9787 (IGizon et alJl2009l iMoradi et al.ll2010h . In order to 
decrease noise, we further average these travel-time maps over a 
region behind the sunspot where scattering phase shifts are large. 
This elliptical portion of the travel-time maps is shown in Fig. [2] 
We use angle brackets (•) to denote this spatial averaging of the 
travel times. 

For one day of observations, the spatially averaged travel- 
time shift introduced by the sunspot is (Aqst) 75 s for the pi 
modes. In order to estimate the noise in these measurements, we 
compute (Aqst) for each of the seven days of observations. We 
define the observed noise level as the standard error in the mean 



5.1. Is it possible to distinguish between sunspot models? 

One way to measure the effect of a change in the sunspot model 
on the waves is to look at the change in the travel-time shifts 
measured with respect to the waves in the quiet-Sun reference 
model, i.e. by computing the difference (Aqst(RS + pert)) - 
<A qs t(RS)). ' 

Table[T]shows these changes in travel-time shifts for the per- 
turbations to the reference sunspot model described in Sects. [23] 
I2.4l and l231 For the case where the Wilson depression is moved 
deeper (shallower) by 50 km, we find that the f wave packet trav- 
els slower (faster) and the pi and P2 wave packets travel faster 
(slower). Because the travel-time shifts for these two cases do 
not have the same amplitude, perturbation theory that is linear in 
the height of the Wilson depression will not fully capture this 
result. The physical cause of these effects is not clear. We 
note, however, that a change in the Wilson depression implies 
a change in the density, sound-speed and Alfv en speed which 
all affect the wave speed dHanasoge et al.ll2012l) . For the case of 
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Table 1. The spatially averaged travel-time shift introduced by the per- 
turbation to the sunspot model, (Aqst(RS + pert)) - (Aqst(RS)). 



Sunspot Perturbation 


Travel-time Shift (s) 

f Pi P2 


Wilson depression: 








50 km deeper 


+0.4 


-2.2 


-3.8 


50 km shallower 


-0.3 


+ 1.9 


+3.2 


Alternative magnetic field 


-2.4 


-2.8 


-3.3 


Sound-speed perturbations: 








A = 0.1, h = -2 Mm 


0.0 


-0.4 


-0.8 


A = 0.1, h = -3 Mm 


0.0 


-0.8 


-0.6 


A = 0.1, h = -4 Mm 


0.0 


-0.8 


-0.2 


A = 0.5,/i = -1.5 Mm 




-1.2 





Notes. Negative (positive) indicates that the waves travel faster (slower) 
through the perturbed model than the reference sunspot model. The 
observational noise level is given in the bottom row. The systematic 
error introduced by the measurement procedure is a 0.1 s. The noise 
level for seven days of SOHO/MDI observations for the f, pi and p2 
mode is 3.7 s, 3.5 s and 1.6 s respectively. The travel-time shifts in bold 
lie above the noise level. 



the p2 wave packet the change in the travel-time shift is roughly 
twice the noise level, implying that the p2 travel times are useful 
for constraining the height of the Wilson depression to an accu- 
racy better than 50 km (assuming the other sunspot parameters 
are well constrained). 

Reducing the strength of the subsurface magnetic field 
causes a change in the travel-time shift of roughly 3 s for all 
wave packets. The physical process governing the interaction of 
the wave with the subsurface magnetic field is not clear. Recall 
that this change in the magnetic field is only below z = -1 Mm. 
As for the case of changes to the height of the Wilson depression, 
for the pi wave packet the change in the travel-time shift is about 
twice that of the observed noise level. Assuming that the Wilson 
depression was well constrained (and all other sunspot parame- 
ters), the subsurface magnetic field may also be constrained. 

Inserting subsurface sound-speed enhancements causes the 
pi and p2 mode wave packets to travel faster than expected. The 
f wave packet is not sensitive to the subsurface sound-speed at 
all, which is also not surprising since it is a surface gravity wave. 
Unless the observed signal to noise can be increased for the p i 
and p2 wave packets by a factor of at least 4, sound-speed per- 
turbations of this type will not be able to be detected with the 
spatially averaged travel-time shifts of these modes. 

5.2. Can we avoid using numerical simulations? 

In the second experiment, we further investigate the sensitivity 
of travel- t ime shifts to changes in th e s ubsurface sound-sp eed. 
iFan et all (I 1 9951) . iKosovichevI i 1 9961) and lBirch et all <2004l) ex- 
plored this problem using perturbation theory around quiet-Sun 
reference models. Here we go beyond small amplitude pertur- 
bation theory and look at finite amplitude perturbations to the 
subsurface sound-speed. In addition, we look at the effects of 
identical sound-speed changes to a reference sunspot model and 
in the quiet-Sun model. 

We insert Gaussian sound-speed perturbations of various 
amplitudes into the reference sunspot model (i.e. Eq. [3] with 
A = -0.01 to 0.5 and h = —1.5 Mm). In order to assess the sig- 
nificance of the reference model, we also insert identical sound- 
speed perturbations into the background quiet-Sun model so that 



co(w, z) = co(z) + Ac(m, z), (5) 

where cq is the quiet-Sun sound speed, and Ac(cj, z) is the pertur- 
bation to the reference sunspot model defined by Eq. [3] We sim- 
ulated the propagation of pi wave packets through these models. 

Figure [3] shows the change in the spatially averaged travel- 
time shifts, (Aqst(RS + 6c)) - (Aqst(RS)), as a function of 
the amplitude of the sound-speed perturbation. The ensuing 
changes to the travel-time shift are not linear with the amplitude 
of the sound-speed perturbation. The deviation from linearity 
for a sound-speed perturbation with a 10% amplitude is 5% - 
a parabolic fit describes the change in the travel-time shifts with 
respect to sound-speed perturbation amplitude better. Notice that 
for all amplitudes, even up to 50%, the change in the travel-time 
shift is well below the observed noise level. Thus, these subsur- 
face sound-speed perturbations are not detectable using only the 
pi wave packet. 

Another way to measure the travel-time shift due to the 
sound-speed perturbations is to directly compute the travel-time 
shift with respect to the reference sunspot model, i.e. (Arst(RS + 
6c)) (see Fig. [3}. This gives a different, but qualitatively similar, 
result to the previous method. This shows that the general be- 
haviour of the change in the travel-time shifts with amplitude is 
not just a consequence of the measurement procedure. 

Figure[3]also shows the travel-time shifts, (Aqst(QS + 6c)), 
caused by imposing the same sound-speed perturbations in the 
quiet-Sun model. In this case, the linear approximation is incor- 
rect by more than 10% when the amplitude of the sound-speed 
perturbation is 10%. More importantly, these travel-time shifts 
are roughly four times larger than the travel-time shifts caused 
by imposing the same sound-speed perturbation in the reference 
sunspot model. As a result, the quiet-Sun model cannot be used 
to predict the travel-time shifts caused by sound-speed pertur- 
bations to a sunspot model because the sunspot substantially 
changes the pi wave packet, and thus its sensitivity to subsur- 
face perturbations. Therefore, it is necessary to use numerical 
simulations for a better understanding of the effects of a sunspot 
on the waves. 



6. Summary and prospects 

We have shown that it should be possible to constrain the height 
of the Wilson depression to within 50 km and to distinguish be- 
tween the two models for the subsurface magnetic field structure 
(a monolithic and a weaker spreading field) using seven days of 
SOHO/MDI observations of one sunspot. Using only the spa- 
tial average of the travel-times shifts it is not possible to detect 
10% sound-speed changes with a length scale of 5 Mm at depths 
between 1.5 and 4 Mm with the same observations. We ex- 
pect that including more information, e.g. the spatial distribution 
of the travel-time shifts, the amplitude of the cross-correlations, 
and higher order modes, will improve the ability of helioseismic 
measurements to constrain the subsurface structure of sunspots. 

We have also demonstrated that the wave travel-time shifts 
caused by changes in sound speed to our sunspot model are dif- 
ferent than the wave travel-time shifts caused by the same change 
in sound speed of the quiet-Sun model. This is due to the ef- 
fect of the sunspot on the wavefield, and this effect is not cap- 
tured by the traditional approximations of local helio seismology 
dFan et al.ll995blKosovichevll996tlBirch et al.l2004l) . Computa- 
tional helioseismology, however, provides a clear path forward. 
These simulations may also be used to test the magnetoghydro- 



Article number, page 4 of [5] 



H. Schunker et al,; Helioseismology of sunspots 




10 20 30 40 50 



A (%) 

Fig. 3. Travel-time perturbations versus sound-speed perturbation am- 
plitudes A, for p[ wave packets. The results from two SLiM numerical 
experiments are shown. In the first procedure, the red diamonds show 
the travel-time perturbations resulting from sound-speed perturbations 
to a background model that includes the reference sunspot, (Aqst(RS + 
6c)) - <Aq S t(RS)). In the second procedure, the black diamonds show 
the travel-time perturbations resulting from the same sound-speed per- 
turbations to a quiet-sun background model, (Aqst(QS + 6c)). In both 
experiments, the travel times (diamonds) are measured using sliding 
quiet-Sun waveforms. For the first experiment, the red asterisks show 
the results using a different measurement procedure, (A RS r(RS + 6c)). 
The solid lines show linear fits to the three smallest perturbations with 
A=-l%, 0, and 1%. The inset shows a zoom around the origin of the 
plot. The dashed curves are parabolic fits. The pi travel-time noise level 
for a 7 days of observation is 3.5 s isindicated by the arrow. 
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Other avenues to explore are the use of vector magnetograms 
and intensity images to further constrain the surface proper- 
ties of sunspots and the use of SDO-HMI Doppler observations 
( I Scherrer et al .1120 121) for helioseismic measurements. 
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